/**
 * @file Mesh.h
 * @author Lishijie (lsj1018845759@outlook.com)
 * @brief 
 * @version 0.1
 * @date 2020-11-08
 * 
 * @copyright Copyright (c) 2020
 * 
 */
#ifndef _MESH_H_
#define _MESH_H_
#include <algorithm>
#include <string>
#include <fstream>
#include <sstream>
#include "Domain.h"
#include "Dofs.h"
//#define unstructured_mesh
#define TEMPLATE template<Dimension DIM>
/**
 * @brief 有限元计算网格类，所有网格从这个类派生；
 *        功能：设置和返回网格大小，计算并返回指定全局编号对应的自由度信息，
 *             返回指定单元的所有自由度的全局编号，以及其他网格信息；
 */
TEMPLATE 
class Mesh
{
public:
    /**
     * @brief 该网格中单元个数；
     * 
     * @return long int 
     */
    virtual long int n_element() = 0;
    /**
     * @brief 该网格中节点个数;
     * 
     * @return long int 
     */
    virtual long int n_mesh_point() = 0;
    /**
     * @brief 该网格中所有自由度总数;
     * 
     * @return long int 
     */
    virtual long int n_dofs() = 0; 
    /**
     * @brief x坐标轴上划分的段数；
     * 
     * @return long& 
     */
    virtual long int& x_seg() = 0;
    /**
     * @brief y坐标轴上划分的段数；
     * 
     * @return long& 
     */
    virtual long int& y_seg() = 0;
    /**
     * @brief z坐标轴上划分的段数；
     * 
     * @return long& 
     */
    virtual long int& z_seg() = 0;
    /**
     * @brief x坐标轴上每一段长度；
     * 
     * @return long& 
     */
    virtual Real x_h() = 0;
    /**
     * @brief y坐标轴上每一段长度；
     * 
     * @return long& 
     */
    virtual Real y_h() = 0;
    /**
     * @brief z坐标轴上每一段长度；
     * 
     * @return long& 
     */
    virtual Real z_h() = 0;
    /** 
     * @brief 每行的自由度个数
     * 
     * @return 
     */
    virtual int n_dofs_per_row() = 0;
    /** 
     * @brief 每列的自由度个数
     * 
     * @return 
     */
    virtual int n_dofs_per_col() = 0;
    /**
     * @brief 返回全局编号为i的自由度
     * 
     * @param i 指定全局编号
     * @return Dofs<DIM> 对应的自由度；
     */
    virtual Dofs<DIM> DofsofIndex(int i) = 0;
    /**
     * @brief 返回指定单元所有全局编号；
     * 
     * @param i 指定单元序数；
     * @return std::vector<int> 存放该单元所有全局编号的数组
     */
    virtual std::vector<int> NodeofEle(int i) = 0;
    /**
     * @brief 返回所有边界上点的全局编号；
     * 
     * @return std::vector<int>& 存放该网格所有边界上自由度全局编号的数组
     */
    virtual std::vector<int>& Boundary() = 0;
    /**
     * @brief 返回计算区域的指针
     * 
     * @return RectangleDomain* 计算区域的指针
     */
    virtual RectangleDomain* RDomain() = 0;
    /**
     * @brief 二维度情况下，第i,j号自由度对应的全局编号
     * 
     * @param i 下标１
     * @param j 下标２
     * @return int 全局编号
     */
    virtual int CortoIdx(int i, int j) = 0;
    /**
     * @brief 判断该index序号点是否一个边界点； 
     * 
     * @param index 需要判断点的全局编号
     * @return true 是边界点
     * @return false 不是边界点
     */
    virtual bool IsBoundary(int index) = 0;
    
    virtual long int n_bndedge(){return 0;};

    virtual std::vector<int> BndEdge(int i){return std::vector<int>(2,0);};
    #ifdef unstructured_mesh
    /**
     * @brief 对于不利用多重网格求解器时进行一些初始化操作；
     */
    virtual void initialize(){};
    #endif
    std::vector<std::vector<int> > _BndEdge;
protected:
    #ifdef unstructured_mesh
    std::vector<Dofs<DIM> > _P;
    std::vector<std::vector<int> > _T;
    std::vector<Dofs<DIM> > _Pb;
    std::vector<std::vector<int> > _Tb;
    std::vector<std::vector<int> > _B;
    std::vector<std::vector<int> > _BndEdge;
    #endif
};
#undef TEMPLATE
/**
 * @brief Q1网格，对应单元为Q1单元；
 * 
 */
/**********************************************************/
/*********************** Q1 2D ****************************/
/**********************************************************/

class Q1Mesh:public Mesh<2>
{
    typedef Mesh<2> _base;
    #ifdef unstructured_mesh
    typedef Mesh<2> _base;
    using _base::_P;
    using _base::_T;
    using _base::_Pb;
    using _base::_Tb;
    using _base::_B;
    #endif
public:
    Q1Mesh(){};
    Q1Mesh(RectangleDomain*);
    Q1Mesh(RectangleDomain*,long int[]);
    Q1Mesh(RectangleDomain*,std::vector<long int>);
    Q1Mesh(std::shared_ptr<RectangleDomain>);
    Q1Mesh(std::shared_ptr<RectangleDomain>,long int[]);
    Q1Mesh(std::shared_ptr<RectangleDomain>,std::vector<long int>);
    ~Q1Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_mesh_point();
    virtual long int n_dofs(); 
    long int n_bndedge(){return (2 * x_seg() + 2 * y_seg());};
    virtual long int& x_seg(){return _seg[0];};
    virtual long int& y_seg(){return _seg[1];}; 
    virtual Real x_h(){return fabs(_RDomain->LeftBottomPoint()[0] - _RDomain->RightBottomPoint()[0]) /static_cast<Real>(x_seg());};
    virtual Real y_h(){return fabs(_RDomain->LeftBottomPoint()[1] - _RDomain->LeftTopPoint()[1]) /static_cast<Real>(y_seg());};
    virtual int n_dofs_per_row(){return _seg[0] + 1;};
    virtual int n_dofs_per_col(){return _seg[1] + 1;};
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
    using _base::_BndEdge;
    #ifdef unstructured_mesh
    virtual void initialize();
    #endif
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;
    virtual long int& z_seg(){return _seg[1];};//二维情况不用，接口就私有吧
    virtual Real z_h(){return y_h();};//二维情况用不到，接口就私有吧
    void ReBoundary();
};
#ifdef unstructured_mesh
void Q1Mesh::initialize()
{
    double xh = x_h();
    double yh = y_h();
    for(int j =  0;j <= y_seg();j++)
        for(int i = 0;i <= x_seg();i++)
        {
            Dofs<2> temp({i * xh,j * yh});
            temp.SetGlobalIndex(CortoIdx(j,i));
            _P.push_back(temp);
        }
    _Pb = _P;
    for(int k = 0;k < n_element();k++)
    {
        std::vector<int> Idx;
        int a = k / (x_seg());
        int b = k % (x_seg());
        Idx.push_back(a * (x_seg() + 1) + b);
        Idx.push_back(Idx[0] + 1);
        Idx.push_back(Idx[1] + (x_seg() + 1));
        Idx.push_back(Idx[2] - 1);
        _T.push_back(Idx);
    }
    _Tb = _T;
    for(int i = 0;i <= _seg[0];i++)
    {
        std::vector<int> b1({-1,i});
        std::vector<int> b2({-1,n_dofs() - i  - 1});
        _B.push_back(b1);
        _B.push_back(b2);
    }
    for(int j = 1;j <= _seg[1] - 1;j++)
    {
        std::vector<int> b3({-1,j*(_seg[0] + 1)});
        std::vector<int> b4({-1,(j + 1)*(_seg[0] + 1) - 1});
        _B.push_back(b3);
        _B.push_back(b4);
    }

    for(int k = 0;k < n_element();k++)
    {
        std::vector<int> idx = _Tb[k];
        for(int i =0;i < 3;i++)
        {
            std::vector<int> edge({-1,k,idx[i],idx[i + 1]});
             _BndEdge.push_back(edge);
        }
        std::vector<int> edge2({-1,k,idx[3],idx[0]});
        _BndEdge.push_back(edge2);
    }
    //auto Ifbndedge = []{std::vector<int> t} ->{return (IsBoundary(t[2]) && IsBoundary(t[3]));};
    std::remove_if(_BndEdge.begin(),_BndEdge.end(),[this](std::vector<int> t){return !(IsBoundary(t[2]) && IsBoundary(t[3]));});
}
#endif
Q1Mesh::Q1Mesh(RectangleDomain* Domain)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

Q1Mesh::Q1Mesh(RectangleDomain* Domain,long int seg[])
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg,seg+2,_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

Q1Mesh::Q1Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif

}

Q1Mesh::Q1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

Q1Mesh::Q1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,long int seg[])
{
     _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
     std::copy(seg,seg+2,_seg);
     ReBoundary();
     #ifdef unstructured_mesh
     initialize();
     #endif
}

Q1Mesh::Q1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

long int Q1Mesh::n_element()
{
    return static_cast<long int>(_seg[0] * _seg[1]);
}
long int Q1Mesh::n_mesh_point()
{
    return static_cast<long int>((_seg[0] + 1) *( _seg[1] + 1));
}
long int Q1Mesh::n_dofs()
{
    return static_cast<long int>((_seg[0] + 1) * (_seg[1] + 1));
}

/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> Q1Mesh::DofsofIndex(int index)
{
    #ifndef unstructured_mesh
    Dofs<2> temp_dofs;
    int i = index % (x_seg() + 1);
    int j = index / static_cast<int>(x_seg() + 1);
    Real x_c = _RDomain->LeftBottomPoint()[0] + i * x_h();
    Real y_c = _RDomain->LeftBottomPoint()[1] + j * y_h();
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    #else
    Dofs<2> temp_dofs = _Pb[index];
    #endif
    return temp_dofs;
} 

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> Q1Mesh::NodeofEle(int i)
{
    #ifndef unstructured_mesh
    std::vector<int> Idx;
    int a = i / (x_seg());
    int b = i % (x_seg());
    Idx.push_back(a * (x_seg() + 1) + b);
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[1] + (x_seg() + 1));
    Idx.push_back(Idx[2] - 1);
    #else
    std::vector<int> Idx = _Tb[i];
    #endif
    return Idx;
}
void Q1Mesh::ReBoundary()
{
    for(int i = 0;i <= _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j*(_seg[0] + 1));
        _BndIndex.push_back((j + 1)*(_seg[0] + 1) - 1);
    }
    //生成边界所有点
     for(int k = 0;k < n_element();k++)
    {
        std::vector<int> idx = NodeofEle(k);
        for(int i =0;i < 3;i++)
        {
            std::vector<int> edge({-1,k,idx[i],idx[i + 1]});
             _BndEdge.push_back(edge);
        }
        std::vector<int> edge2({-1,k,idx[3],idx[0]});
        _BndEdge.push_back(edge2);
    }//默认边界条件为Dirichlet;
    std::remove_if(_BndEdge.begin(),_BndEdge.end(),[this](std::vector<int> t){return !(IsBoundary(t[2]) && IsBoundary(t[3]));});
    _BndEdge.resize(2 * x_seg() + 2 * y_seg());
}

/**
 * @brief P1网格，对应单元为P1单元；
 * 
 */
/**********************************************************/
/*********************** P1 2D ****************************/
/**********************************************************/

/* 如何划分单元格？                     
 * 首先将实际计算区域划分为规则的矩形网格，对于每个矩形网格，
 * 将其第0个和第2个顶点连线，划分为两个三角形，从左至右排序。
 * 2  ——  3 
 * |     /|    如图左所示，0-3-2, 0-1-3 构成两个 P1 单元片，
 * |    / |    且 0-3-2 编号在前。
 * |   /  |
 * |  /   |
 * | /    |
 * |/     |
 * 0  ——  1
 *
 */

class P1Mesh:public Mesh<2>
{
public:
    #ifdef unstructured_mesh
    typedef Mesh<2> _base;
    using _base::_P;
    using _base::_T;
    using _base::_Pb;
    using _base::_Tb;
    using _base::_B;
    using _base::_BndEdge;
    #endif
public:
    P1Mesh(){};
    P1Mesh(RectangleDomain*);
    P1Mesh(RectangleDomain*, long int[]);
    P1Mesh(RectangleDomain*, std::vector<long int>);
    P1Mesh(std::shared_ptr<RectangleDomain>);
    P1Mesh(std::shared_ptr<RectangleDomain>, long int[]);
    P1Mesh(std::shared_ptr<RectangleDomain>, std::vector<long int>);
    ~P1Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_mesh_point();
    virtual long int n_dofs(); 
    virtual long int& x_seg(){return _seg[0];};
    virtual long int& y_seg(){return _seg[1];}; 
    virtual Real x_h(){return fabs(_RDomain->LeftBottomPoint()[0] - _RDomain->RightBottomPoint()[0]) / static_cast<Real>(x_seg());};
    virtual Real y_h(){return fabs(_RDomain->LeftBottomPoint()[1] - _RDomain->LeftTopPoint()[1]) / static_cast<Real>(y_seg());};
    virtual int n_dofs_per_row(){return _seg[0] + 1;};
    virtual int n_dofs_per_col(){return _seg[1] + 1;};
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
    #ifdef unstructured_mesh
    virtual void initialize();
    #endif
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;
    virtual long int& z_seg(){return _seg[1];};
    virtual Real z_h(){return y_h();};
    void ReBoundary();
};
#ifdef unstructured_mesh
void P1Mesh::initialize()
{
    double xh = x_h();
    double yh = y_h();
    for(int j =  0;j <= y_seg();j++)
        for(int i = 0;i <= x_seg();i++)
        {
            Dofs<2> temp({i * xh,j * yh});
            temp.SetGlobalIndex(CortoIdx(j,i));
            _P.push_back(temp);
        }
    _Pb = _P;
    for(int i = 0;i < n_element();i++)
    {
        std::vector<int> Idx;
        int a = i / (2 * x_seg());
        int b = i % (2 * x_seg());
        if (b % 2 == 0)
        {
	        Idx.push_back(a * (x_seg() + 1) + b / 2);
	        Idx.push_back(Idx[0] + x_seg() + 2);
	        Idx.push_back(Idx[1] - 1);
            _T.push_back(Idx);
        }
        else
        {
	        Idx.push_back(a * (x_seg() + 1) + b / 2);
	        Idx.push_back(Idx[0] + 1);
	        Idx.push_back(Idx[1] + x_seg() + 1);
            _T.push_back(Idx);
        }
    }
    _Tb = _T;
}
#endif
P1Mesh::P1Mesh(RectangleDomain* Domain)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

P1Mesh::P1Mesh(RectangleDomain* Domain, long int seg[])
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg, seg+2, _seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

P1Mesh::P1Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif

}

P1Mesh::P1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

P1Mesh::P1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,long int seg[])
{
     _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
     std::copy(seg,seg+2,_seg);
     ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

P1Mesh::P1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

long int P1Mesh::n_element()
{
    return static_cast<long int>(2 * _seg[0] * _seg[1]);
}
long int P1Mesh::n_mesh_point()
{
    return static_cast<long int>((_seg[0] + 1) *( _seg[1] + 1));
}
long int P1Mesh::n_dofs()
{
    return static_cast<long int>((_seg[0] + 1) * (_seg[1] + 1)) ;
}

/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> P1Mesh::DofsofIndex(int index)
{
    #ifndef unstructured_mesh
    Dofs<2> temp_dofs;
    int i = index % (x_seg() + 1);
    int j = index / static_cast<int>(x_seg() + 1);
    Real x_c = _RDomain->LeftBottomPoint()[0] + i * x_h();
    Real y_c = _RDomain->LeftBottomPoint()[1] + j * y_h();
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    #else
    Dofs<2> temp_dofs = _Pb[index];
    #endif
    return temp_dofs;
}

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> P1Mesh::NodeofEle(int i)
{
    #ifndef unstructured_mesh
    std::vector<int> Idx;
    int a = i / (2 * x_seg());
    int b = i % (2 * x_seg());
    if (b % 2 == 0)
    {
	Idx.push_back(a * (x_seg() + 1) + b / 2);
	Idx.push_back(Idx[0] + x_seg() + 2);
	Idx.push_back(Idx[1] - 1);
    }
    else
    {
	Idx.push_back(a * (x_seg() + 1) + b / 2);
	Idx.push_back(Idx[0] + 1);
	Idx.push_back(Idx[1] + x_seg() + 1);
    }
    #else
    std::vector<int> Idx = _Tb[i];
    #endif
    return Idx;
}

void P1Mesh::ReBoundary()
{
    for(int i = 0;i <= _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j*(_seg[0] + 1));
        _BndIndex.push_back((j + 1)*(_seg[0] + 1) - 1);
    }
}
/**
 * @brief Q2网格，对应单元为Q2单元；
 * 
 */
/**********************************************************/
/*********************** Q2 2D ****************************/
/**********************************************************/
class Q2Mesh:public Mesh<2>
{
public:
    #ifdef unstructured_mesh
    typedef Mesh<2> _base;
    using _base::_P;
    using _base::_T;
    using _base::_Pb;
    using _base::_Tb;
    using _base::_B;
    using _base::_BndEdge;
    #endif
public:

    Q2Mesh(){};
    Q2Mesh(RectangleDomain*);
    Q2Mesh(RectangleDomain*,long int[]);
    Q2Mesh(RectangleDomain*,std::vector<long int>);
    Q2Mesh(std::shared_ptr<RectangleDomain>);
    Q2Mesh(std::shared_ptr<RectangleDomain>,long int[]);
    Q2Mesh(std::shared_ptr<RectangleDomain>,std::vector<long int>);
    ~Q2Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_mesh_point();
    virtual long int n_dofs(); 
    virtual long int& x_seg(){return _seg[0];};
    virtual long int& y_seg(){return _seg[1];}; 
    virtual Real x_h(){return fabs(_RDomain->LeftBottomPoint()[0] - _RDomain->RightBottomPoint()[0]) /static_cast<Real>(x_seg());};
    virtual Real y_h(){return fabs(_RDomain->LeftBottomPoint()[1] - _RDomain->LeftTopPoint()[1]) /static_cast<Real>(y_seg());};
    virtual int n_dofs_per_row(){return 2*_seg[0] + 1;};
    virtual int n_dofs_per_col(){return 2*_seg[1] + 1;};
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(2*_seg[0]+1) + j);};//Q2 每行有两倍多的点。
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
    #ifdef unstructured_mesh
    virtual void initialize();
    #endif
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;
    virtual long int& z_seg(){return _seg[1];};//二维情况不用，接口就私有吧
    virtual Real z_h(){return y_h();};//二维情况用不到，接口就私有吧
    void ReBoundary();
};

#ifdef unstructured_mesh
void Q2Mesh::initialize()
{
    double xh = x_h();
    double yh = y_h();
    for(int j =  0;j <= y_seg();j++)
        for(int i = 0;i <= x_seg();i++)
        {
            Dofs<2> temp({i * xh,j * yh});
            temp.SetGlobalIndex(j *  2 * (2 * x_seg() + 1) + i * 2);
            _P.push_back(temp);
        }

    for(int j =  0;j <= 2 * y_seg();j++)
        for(int i = 0;i <= 2 * x_seg();i++)
        {
            Dofs<2> tempb({i * xh / 2.0,j * yh / 2.0});
            tempb.SetGlobalIndex(CortoIdx(j,i));
            _Pb.push_back(tempb);
        }
    for(int i = 0;i < n_element();i++)
    {
        std::vector<int> Idx;
        std::vector<int> Idxb;
        int a = i / (x_seg());
        int b = i % (x_seg());
        Idx.push_back(a * 2 * (2 * x_seg() + 1) + 2 * b);
        Idx.push_back(Idx[0] + 2);
        Idx.push_back(Idx[1] + 2 * (2 * x_seg() + 1));
        Idx.push_back(Idx[2] - 2);
        Idxb.push_back(a * 2 * (2 * x_seg() + 1) + 2 * b);
        Idxb.push_back(Idxb[0] + 2);
        Idxb.push_back(Idxb[1] + 2 * (2 * x_seg() + 1));
        Idxb.push_back(Idxb[2] - 2);
        Idxb.push_back(Idxb[0] + 1);
        Idxb.push_back(Idxb[1] + (2 * x_seg() + 1));
        Idxb.push_back(Idxb[2] - 1);
        Idxb.push_back(Idxb[0] + (2 * x_seg() + 1));
        Idxb.push_back(Idxb[7] + 1);
        _T.push_back(Idx);
        _Tb.push_back(Idxb);
    }
    
}
#endif
Q2Mesh::Q2Mesh(RectangleDomain* Domain)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

Q2Mesh::Q2Mesh(RectangleDomain* Domain,long int seg[])
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg,seg+2,_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

Q2Mesh::Q2Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif

}

Q2Mesh::Q2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

Q2Mesh::Q2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,long int seg[])
{
     _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
     std::copy(seg,seg+2,_seg);
     ReBoundary();
     #ifdef unstructured_mesh
     initialize();
     #endif
}

Q2Mesh::Q2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

long int Q2Mesh::n_element()
{
    return static_cast<long int>(_seg[0] * _seg[1]);
}
long int Q2Mesh::n_mesh_point()
{
    return static_cast<long int>((_seg[0] + 1) *( _seg[1] + 1));
}
long int Q2Mesh::n_dofs()
{
    return static_cast<long int>((2 * _seg[0] + 1) * (2 * _seg[1] + 1));
}

/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> Q2Mesh::DofsofIndex(int index)
{
    #ifndef unstructured_mesh
    Dofs<2> temp_dofs;
    int i = index % (2 * x_seg() + 1);
    int j = index / static_cast<int>(2 * x_seg() + 1);
    Real x_c = _RDomain->LeftBottomPoint()[0] + i * x_h() /2.0;
    Real y_c = _RDomain->LeftBottomPoint()[1] + j * y_h() /2.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    #else
    Dofs<2> temp_dofs = _Pb[index];
    #endif
    return temp_dofs;
} 

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> Q2Mesh::NodeofEle(int i)
{
    #ifndef unstructured_mesh
    std::vector<int> Idx;
    int a = i / (x_seg());
    int b = i % (x_seg());
    Idx.push_back(a * 2 * (2 * x_seg() + 1) + 2 * b);
    Idx.push_back(Idx[0] + 2);
    Idx.push_back(Idx[1] + 2 * (2 * x_seg() + 1));
    Idx.push_back(Idx[2] - 2);
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[1] + (2 * x_seg() + 1));
    Idx.push_back(Idx[2] - 1);
    Idx.push_back(Idx[0] + (2 * x_seg() + 1));
    Idx.push_back(Idx[7] + 1);
    #else
    std::vector<int> Idx = _Tb[i];
    #endif
    return Idx;
}
void Q2Mesh::ReBoundary()
{
    for(int i = 0;i <= 2 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= 2 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (2  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(2 * _seg[0] + 1) - 1);
    }
}
/**
 * @brief P2网格，对应单元为P2单元；
 * 
 */
/**********************************************************/
/*********************** P2 2D ****************************/
/**********************************************************/

/* 如何划分单元格？                     
 * 首先将实际计算区域划分为规则的矩形网格，对于每个矩形网格，
 * 将其第0个和第8个顶点连线，划分为两个三角形，从左至右排序。
 * 从左下角为第一个号，逆时针排序．
 * 6  —7—  8 
 * |     /|    如图左所示，0-8-6-4-7-3, 0-2-8-1-5-4 构成两个 P2 单元片，
 * |    / |    且 0-8-6-4-7-3 编号在前。
 * 3   4  5
 * |  /   |
 * | /    |
 * |/     |
 * 0  —1—  2
 *
 */

class P2Mesh:public Mesh<2>
{
public:
    #ifdef unstructured_mesh
    typedef Mesh<2> _base;
    using _base::_P;
    using _base::_T;
    using _base::_Pb;
    using _base::_Tb;
    using _base::_B;
    using _base::_BndEdge;
    #endif
public:
    P2Mesh(){};
    P2Mesh(RectangleDomain*);
    P2Mesh(RectangleDomain*, long int[]);
    P2Mesh(RectangleDomain*, std::vector<long int>);
    P2Mesh(std::shared_ptr<RectangleDomain>);
    P2Mesh(std::shared_ptr<RectangleDomain>, long int[]);
    P2Mesh(std::shared_ptr<RectangleDomain>, std::vector<long int>);
    ~P2Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_mesh_point();
    virtual long int n_dofs(); 
    virtual long int& x_seg(){return _seg[0];};
    virtual long int& y_seg(){return _seg[1];}; 
    virtual Real x_h(){return fabs(_RDomain->LeftBottomPoint()[0] - _RDomain->RightBottomPoint()[0]) / static_cast<Real>(x_seg());};
    virtual Real y_h(){return fabs(_RDomain->LeftBottomPoint()[1] - _RDomain->LeftTopPoint()[1]) / static_cast<Real>(y_seg());};
    virtual int n_dofs_per_row(){return 2*_seg[0] + 1;};
    virtual int n_dofs_per_col(){return 2*_seg[1] + 1;};
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(2*_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
    #ifdef unstructured_mesh
    virtual void initialize();
    #endif
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;
    virtual long int& z_seg(){return _seg[1];};
    virtual Real z_h(){return y_h();};
    void ReBoundary();
};
#ifdef unstructured_mesh
void P2Mesh::initialize()
{
    double xh = x_h();
    double yh = y_h();
    for(int j =  0;j <= y_seg();j++)
        for(int i = 0;i <= x_seg();i++)
        {
            Dofs<2> temp({i * xh,j * yh});
            temp.SetGlobalIndex(j *  2 * (2 * x_seg() + 1) + i * 2);
            _P.push_back(temp);
        }

    for(int j =  0;j <= 2 * y_seg();j++)
        for(int i = 0;i <= 2 * x_seg();i++)
        {
            Dofs<2> tempb({i * xh / 2.0,j * yh / 2.0});
            tempb.SetGlobalIndex(CortoIdx(j,i));
            _Pb.push_back(tempb);
        }
    for(int i = 0;i < n_element();i++)
    {
        std::vector<int> Idx;
        std::vector<int> Idxb;
        int a = i / (2 * x_seg());
        int b = i % (2 * x_seg());
        if (b % 2 == 0)
        {
        Idx.push_back(a * 2 * (2 * x_seg() + 1) + b);
	    Idx.push_back(Idx[0] + 2 * (2 * x_seg() + 1) + 2);
	    Idx.push_back(Idx[1] - 2);
	    Idxb.push_back(a * 2 * (2 * x_seg() + 1) + b);
	    Idxb.push_back(Idxb[0] + 2 * (2 * x_seg() + 1) + 2);
	    Idxb.push_back(Idxb[1] - 2);
        Idxb.push_back(Idxb[0] + 1 + (2 * x_seg() + 1));
        Idxb.push_back(Idxb[1] - 1);
        Idxb.push_back(Idxb[3] - 1);
        _T.push_back(Idx);
        _Tb.push_back(Idxb);
        }
        else
        {
        Idx.push_back(a * 2 *(2 * x_seg() + 1) + b - 1);
	    Idx.push_back(Idx[0] + 2);
	    Idx.push_back(Idx[1] + 2 * (2 * x_seg() + 1));
	    Idxb.push_back(a * 2 *(2 * x_seg() + 1) + b - 1);
	    Idxb.push_back(Idxb[0] + 2);
	    Idxb.push_back(Idxb[1] + 2 * (2 * x_seg() + 1));
        Idxb.push_back(Idxb[0] + 1);
        Idxb.push_back(Idxb[1] + 2 *x_seg() + 1);
        Idxb.push_back(Idxb[4] - 1);
        _T.push_back(Idx);
        _Tb.push_back(Idxb);
        }
    }
}
#endif
P2Mesh::P2Mesh(RectangleDomain* Domain)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

P2Mesh::P2Mesh(RectangleDomain* Domain, long int seg[])
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg, seg+2, _seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

P2Mesh::P2Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif

}

P2Mesh::P2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

P2Mesh::P2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,long int seg[])
{
     _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
     std::copy(seg,seg+2,_seg);
     ReBoundary();
     #ifdef unstructured_mesh
    initialize();
    #endif
}

P2Mesh::P2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
    #ifdef unstructured_mesh
    initialize();
    #endif
}

long int P2Mesh::n_element()
{
    return static_cast<long int>(2 * _seg[0] * _seg[1]);
}
long int P2Mesh::n_mesh_point()
{
    return static_cast<long int>((_seg[0] + 1) *( _seg[1] + 1));
}
long int P2Mesh::n_dofs()
{
    return static_cast<long int>((2 * _seg[0] + 1) * (2 * _seg[1] + 1)) ;
}

/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> P2Mesh::DofsofIndex(int index)
{
    #ifndef unstructured_mesh
    Dofs<2> temp_dofs;
    int i = index % (2 * x_seg() + 1);
    int j = index / static_cast<int>(2 * x_seg() + 1);
    Real x_c = _RDomain->LeftBottomPoint()[0] + i * x_h() / 2.0;
    Real y_c = _RDomain->LeftBottomPoint()[1] + j * y_h() / 2.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    #else
    Dofs<2> temp_dofs = _Pb[index];
    #endif
    return temp_dofs;
}

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> P2Mesh::NodeofEle(int i)
{
    #ifndef unstructured_mesh
    std::vector<int> Idx;
    int a = i / (2 * x_seg());
    int b = i % (2 * x_seg());
    if (b % 2 == 0)
    {
	Idx.push_back(a * 2 * (2 * x_seg() + 1) + b);
	Idx.push_back(Idx[0] + 2 * (2 * x_seg() + 1) + 2);
	Idx.push_back(Idx[1] - 2);
    Idx.push_back(Idx[0] + 1 + (2 * x_seg() + 1));
    Idx.push_back(Idx[1] - 1);
    Idx.push_back(Idx[3] - 1);
    }
    else
    {
	Idx.push_back(a * 2 *(2 * x_seg() + 1) + b - 1);
	Idx.push_back(Idx[0] + 2);
	Idx.push_back(Idx[1] + 2 * (2 * x_seg() + 1));
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[1] + 2 *x_seg() + 1);
    Idx.push_back(Idx[4] - 1);
    }
    #else
    std::vector<int> Idx = _Tb[i];
    #endif
    return Idx;
}

void P2Mesh::ReBoundary()
{
    for(int i = 0;i <= 2 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_dofs()- i - 1);
    }
    for(int j = 1;j <= 2 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (2  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(2 * _seg[0] + 1) - 1);
    }
}

/**********************************************************/
/********************* Eazymesh 2D ************************/
/**********************************************************/

/// 基于 Eazymesh 数据建立的派生类 EazyMesh。
class EazyMesh:public Mesh<2>
{
    typedef std::string string;
    typedef std::stringstream stringstream;
    typedef std::ifstream ifstream;
private:
    string _MeshName;
    long int _ndofs;
    long int _nele;
    long int _nbndedge;
    std::vector<Dofs<2> > _DofsList;
    std::vector<int> _BndDofs;
    std::vector<std::vector<int> > _ElementData;
    std::vector<std::vector<int> > _BndEdges; 

    /// 暂时不用，私有
    long int a;
    virtual long int& x_seg(){return a;};
    virtual long int& y_seg(){return a;}; 
    virtual Real x_h(){return -1;};
    virtual Real y_h(){return -1;};
    virtual int n_dofs_per_row(){return -1;};
    virtual int n_dofs_per_col(){return -1;};
    virtual RectangleDomain* RDomain() {return 0;};
    virtual int CortoIdx(int i, int j) {return -1;};
    virtual long int& z_seg(){return a;};
    virtual Real z_h(){return -1;};
    
public:
    EazyMesh();
    EazyMesh(string meshname);
    void InputDofData();
    void InputEdgeData();
    void InputEleData();
    virtual long int n_element(){return _nele;};
    virtual long int n_mesh_point(){return _ndofs;};
    virtual long int n_dofs(){return _ndofs;};
    long int n_bndedge(){return _nbndedge;};
    virtual Dofs<2> DofsofIndex(int i){return _DofsList[i];};
    virtual std::vector<int> NodeofEle(int i){return _ElementData[i];};
    virtual std::vector<int>& Boundary(){return _BndDofs;};
    virtual bool IsBoundary(int index){return !(std::find(_BndDofs.begin(),_BndDofs.end(),index)==_BndDofs.end());};
    std::vector<int> BndEdge(int i){return _BndEdges[i];};
    

    
};

EazyMesh::EazyMesh(string meshname)
{
    _MeshName = meshname;
    InputDofData();
    InputEdgeData();
    InputEleData();
}

void EazyMesh::InputDofData()
{
    ifstream fin;
    fin.open(_MeshName+".n");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    ss >> _ndofs;
    _DofsList.resize(_ndofs);

    int a, bndflag;
    double x, y;
    for (int i = 0; i < _ndofs; i++)
    {
	getline(fin, line);
	ss.clear();
	ss.str(line);
	ss >> a >> x >> y >> bndflag;
	Dofs<2> temp_dof;
	temp_dof[0] = x;
	temp_dof[1] = y;
	temp_dof.GlobalIndex() = i;
	_DofsList[i] = temp_dof;
	if (bndflag)
	    _BndDofs.push_back(i);	    
    }
}

void EazyMesh::InputEdgeData()
{
    ifstream fin;
    fin.open(_MeshName+".s");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    int _nedge;
    ss >> _nedge;
    int a, d0, d1, bndflag;
    for (int i = 0; i < _nedge; i++)
    {
	getline(fin, line);
	ss.clear();
	ss.str(line);
	ss >> a >> d0 >> d1 >> a >> a >> bndflag;
	if (bndflag)
	{
	    std::vector<int> temp_edge(2);
	    temp_edge[0] = d0;
	    temp_edge[1] = d1;
	    _BndEdges.push_back(temp_edge);
 	}
    }
    _nbndedge = _BndEdges.size();
}

void EazyMesh::InputEleData()
{
    ifstream fin;
    fin.open(_MeshName+".e");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    ss >> _nele;
    _ElementData.resize(_nele);

    int a, d0, d1 ,d2;
    for (int i = 0; i < _nele; i++)
    {
	getline(fin, line);
	ss.clear();
	ss.str(line);
	ss >> a >> d0 >> d1 >> d2;
	_ElementData[i].resize(3);
	_ElementData[i][0] = d0;
	_ElementData[i][1] = d1;
	_ElementData[i][2] = d2;
    }
}

#else
#endif
